
clc
clear all
tic

%%%
%%% Step 1. Calibration under f=0
%%%

% variable set-up

params.alpha=0.5;
params.eta=6;   
params.r=0.02;
params.pi=0.001;

guess_cal=[0.024 0.438 0.114 0.092 459.667 161.327 817.894];

params.weight=[1 1 1 1 1 1 1 1 1];
params.weight=params.weight./sum(params.weight);

iter=0;
dist=1;
while dist>0.001 % average percent deviation 0.1%
    iter=iter+1
% local min
    [xx_cal,fval_cal] = fminsearchbnd(@(x) calibration_function_f0(x,params),guess_cal,[0 0 0 0 0 0 0],[1 1 1 1 10^4 10^4 10^4]);
    dist=0;
    for i=1:7
        dist=dist+(1/7)*(abs(xx_cal(i)-guess_cal(i))/abs(guess_cal(i)));
    end
    dist
    guess_cal=xx_cal;
end

% simulation
    params.lambdaE=xx_cal(1);
    params.lambdaS=xx_cal(2);
    params.lambdaI=xx_cal(3);
    params.lambdaC=xx_cal(4);
    params.f=0;
    params.xiE=xx_cal(5);
    params.xiI=xx_cal(6);
    params.lambdae=params.lambdaE;
    params.xie=xx_cal(7);

    guess_sim=0.1*ones(10,1);
    options = optimset('MaxFunEvals',10^5,'MaxIter',10^5,'TolFun',1e-30);
    [xx_sim,fval_sim] = fminsearchbnd(@(x) equation_difference(x,params),guess_sim,[0 0 0 0 0 0 0 0 0 0],[],options);
    
    alpha=params.alpha;
    pi=params.pi;
    eta=params.eta;
    r=params.r;
    lambdaE=params.lambdaE;
    lambdaS=params.lambdaS;
    lambdaI=params.lambdaI;
    lambdaC=params.lambdaC;
    lambdae=params.lambdae;
    f=params.f;
    xiE=params.xiE;
    xiI=params.xiI;
    xie=params.xie;
    
    xI=xx_sim(2);
    xE=xx_sim(5);
    xe=xx_sim(8);
    tau=xE+xe;
    theta=xe/(xe+xE+xI);
    B=xx_sim(6);
    gq=xx_sim(7);
    RD_all=xiI*eta^(-1)*(xI)^(1/(1-alpha))+xiE*eta^(-1)*(xE)^(1/(1-alpha))+xie*eta^(-1)*(xe)^(1/(1-alpha))+B*xI*eta^(-1)*pi^(-1)*(1-f);
    RD_new=xie*eta^(-1)*(xe)^(1/(1-alpha));
    
    moment1=xI*(eta-1)*log(1-lambdaC)+(xE+xe)*(eta-1)*log(1-lambdaS); % target -0.1506
    moment2=((xE+xe)*log(1-lambdaS))/(xI*log(1-lambdaC)); % target 1.4858
    moment3=xE*((1+lambdaE)^(eta-1)+(1-lambdaS)^(eta-1)-2)+xI*((1+lambdaI)^(eta-1)+(1-lambdaC)^(eta-1)-2)+xe*((1+lambdae)^(eta-1)+(1-lambdaS)^(eta-1)-2); % target 0.016
    moment4=xe*(1+lambdae)^(eta-1); % target 0.0114
    moment5=xI/xE; % target 6
    moment6=RD_all; % target 0.0988
    moment7=xI+xE; % target 0.145
    moment8=xe; % target 0.0108
    moment9=RD_new; % target 0.0115
    
    moments=[moment1 moment2 moment3 moment4 moment5 moment6 moment7 moment8 moment9]
    
    weight=params.weight;
    yyy=weight(1)*(abs(moment1+0.1506)/0.1506) + weight(2)*(abs(moment2-1.4858)/1.4858) + weight(3)*(abs(moment3-0.016)/0.016) + weight(4)*(abs(moment4-0.0114)/0.0114) + weight(6)*(abs(moment6-0.0988)/0.0988) + weight(7)*(abs(moment7-0.145)/0.145) + weight(8)*(abs(moment8-0.0108)/0.0108) + weight(9)*(abs(moment9-0.0115)/0.0115);

    fval_cal
    fval_sim
    
    xx_eqm=[xI xE xe gq]
   
    save calibration_outcome_f0.mat xx_cal xx_eqm
    iter

   
%%%
%%% Step 2. Comparing Two Cases (f>0, f=0)
%%%

clear all

% variable set-up

params.alpha=0.5;
params.eta=6;  
params.r=0.02;
params.pi=0.001;

%% Load Calibration Outcome
load calibration_outcome.mat

f_mat=[xx_cal(5); 0];
S=size(f_mat,1);
L=80;
xiE_mat=linspace(50,3000,L);
xI_mat=zeros(L,S);

params.lambdaE=xx_cal(1);
params.lambdaS=xx_cal(2);

params.lambdaI=xx_cal(3);
params.lambdaC=xx_cal(4);

params.lambdae=params.lambdaE;

params.xiI=xx_cal(7);
params.xie=xx_cal(8);

for s=1:1
for l=1:L
    params.f=f_mat(s);
    params.xiE=xiE_mat(l);
    
    guess_sim=0.1*ones(10,1);
    options = optimset('MaxFunEvals',10^5,'MaxIter',10^5,'TolFun',1e-30);
    [xx_sim,fval_sim] = fminsearchbnd(@(x) equation_difference(x,params),guess_sim,[0 0 0 0 0 0 0 0 0 0],[],options);
    xx_sim;
    fval_sim
    xI=xx_sim(2);
    xE=xx_sim(5);
    
    if fval_sim<1e-15
    xI_mat(l,s)=xI;
    else
        xI_mat(l,s)=NaN;
    end
end
end

%% Load Calibration Outcome
load calibration_outcome_f0.mat

params.lambdaE=xx_cal(1);
params.lambdaS=xx_cal(2);

params.lambdaI=xx_cal(3);
params.lambdaC=xx_cal(4);

params.lambdae=params.lambdaE;

params.xiI=xx_cal(6);
params.xie=xx_cal(7);


for s=2:2
for l=1:L
    params.f=f_mat(s);
    params.xiE=xiE_mat(l);
    
    guess_sim=0.1*ones(10,1);
    options = optimset('MaxFunEvals',10^5,'MaxIter',10^5,'TolFun',1e-30);
    [xx_sim,fval_sim] = fminsearchbnd(@(x) equation_difference(x,params),guess_sim,[0 0 0 0 0 0 0 0 0 0],[],options);
    xx_sim;
    fval_sim
    xI=xx_sim(2);
    xE=xx_sim(5);
    
    if fval_sim<1e-15
    xI_mat(l,s)=xI;
    else
        xI_mat(l,s)=NaN;
    end
end
end

figure(1)
hold on
plot(xiE_mat,xI_mat(:,1),'LineWidth',4,'Color','black')
plot(xiE_mat,xI_mat(:,2),'--','LineWidth',4,'Color',[.5 .5 .5])
hold off
xline(459.666792678737,'-',{'Calibrated','f=17.6'});
xline(2402.77657186523,'-',{'Calibrated','f=0'});
legend({'f=17.6','f=0'},'FontSize',14,'Location','north')
xlabel('External Shifter Parameter');
ylabel('Likelihood of Internal Product Improvement');
grid on;
exportgraphics(gcf,'OUTPUT/Figure9.eps','BackgroundColor','none','ContentType','vector')
